Homework 5
Describe the setup and each step in your solutions with words and clearly label your final answers. Use Matlab for plotting and programming and include your code as an appendix to your problem set.
Table of Contents
Problem 1
Streeter–Phelps Model. The Streeter–Phelps equation is used in the study of water pollution as a water quality modeling tool. The model describes how dissolved oxygen (DO) decreases in a river or stream along a certain distance by degradation of biochemical oxygen demand (BOD). The equation was derived by H. W. Streeter, a sanitary engineer, and Earle B. Phelps, a consultant for the U.S. Public Health Service, in 1925, based on field data from the Ohio River. The equation is also known as the DO sag equation.
In this problem set, we consider the Streeter-Phelps model in a 2D river. The governing equations for BOD and DO includes advection and diffusion process, with additional source and sink terms.
∂ B ∂ t + u i ∂ B ∂ x i = D i ∂ 2 B ∂ x i 2 − K r B + δ ( x − x 0 , y − y 0 ) Q B , \begin{equation*}
\tag{1}
\frac{\partial B}{\partial t} + u_i \frac{\partial B}{\partial x_i} = D_i \frac{\partial^2 B}{\partial x_i^2} - K_{\text{r}} B + \delta(x - x_0, y - y_0) Q_{\text{B}},
\end{equation*} ∂ t ∂ B + u i ∂ x i ∂ B = D i ∂ x i 2 ∂ 2 B − K r B + δ ( x − x 0 , y − y 0 ) Q B , ( 1 )
∂ O ∂ t + u i ∂ O ∂ x i = D i ∂ 2 O ∂ x i 2 + K a ( O sat − O ) − K r B , \begin{equation*}
\tag{2}
\frac{\partial O}{\partial t} + u_i \frac{\partial O}{\partial x_i} = D_i \frac{\partial^2 O}{\partial x_i^2} + K_{\text{a}} (O_{\text{sat}} - O) - K_{\text{r}} B,
\end{equation*} ∂ t ∂ O + u i ∂ x i ∂ O = D i ∂ x i 2 ∂ 2 O + K a ( O sat − O ) − K r B , ( 2 )
where B B B and O O O are BOD and DO concentrations with units [ M L − 3 ] [\mathrm{M} \mathrm{L}^{−3}] [ M L − 3 ] . K r K_{\text{r}} K r [ T − 1 ] [\mathrm{T}^{−1}] [ T − 1 ] is the reaction rate between BOD and DO. K a K_{\text{a}} K a [ T − 1 ] [\mathrm{T}^{−1}] [ T − 1 ] is the mass transfer coefficient for oxygen. The reaeration rate is given by K a ( O sat − O ) K_{\text{a}} (O_{\text{sat}} −O) K a ( O sat − O ) , where O sat O_{\text{sat}} O sat is the saturation oxygen content in water exposed to air, it is mostly affected by water temperature. The reaeration term describes the flux of oxygen across the air-water interface into the water. Q B Q_{\text{B}} Q B [ M T − 1 ] [\mathrm{M} \mathrm{T}^{−1}] [ M T − 1 ] is the rate of BOD discharge into the river. δ ( x − x 0 , y − y 0 ) \delta(x −x_0, y −y_0) δ ( x − x 0 , y − y 0 ) is a delta function indicating a point source. i = 1 , 2 i = 1, 2 i = 1 , 2 for the x x x and y y y directions. Repeated index suggest summation over each value of the repeated index, e.g.,
u i ∂ B ∂ x i ≡ u 1 ∂ B ∂ x 1 + u 2 ∂ B ∂ x 2 = u ∂ B ∂ x + v ∂ B ∂ y . \begin{equation*}
u_i \frac{\partial B}{\partial x_i} \equiv u_1 \frac{\partial B}{\partial x_1} + u_2 \frac{\partial B}{\partial x_2} = u \frac{\partial B}{\partial x} + v \frac{\partial B}{\partial y}.
\end{equation*} u i ∂ x i ∂ B ≡ u 1 ∂ x 1 ∂ B + u 2 ∂ x 2 ∂ B = u ∂ x ∂ B + v ∂ y ∂ B .
1. Discretize the above two equations using the forward in time, centered in space (FTCS) scheme. Set the proper lateral boundary conditions for the problem (HINT: use solid wall boundaries in the cross flow y y y direction, set Neumann (zero-gradient) conditions for the outflow boundary, and Dirichlet conditions for the inflow boundary. Try using the ghost point technique for setting boundary conditions).
2. Perform von Neumann stability analysis for the discretized equations. Determine the proper proper Courant number Cr = u Δ t / Δ x , r = D Δ t / Δ x 2 \text{Cr} = u \Delta t / \Delta x, \, r = D \Delta t / {\Delta x}^2 Cr = u Δ t /Δ x , r = D Δ t / Δ x 2 , and Peclet number Pe = u Δ x / D \text{Pe} = u \Delta x / D Pe = u Δ x / D values for a stable solution (HINT: you may combine the two equations, by ignoring the BOD discharge term δ ( x − x 0 , y − y 0 ) Q B \delta(x −x_0, y − y_0) Q_{\text{B}} δ ( x − x 0 , y − y 0 ) Q B in Eq. (1), and the K a O sat − K r B K_{\text{a}} O_{\text{sat}} − K_{\text{r}} B K a O sat − K r B terms in Eq. (2). You may further simply the analysis by assuming D x = D y D_x = D_y D x = D y ).
3. Perform modified equation analysis, and discuss the truncation errors.
4. Write a Matlab script to solve the problem. Use the following assumptions:
The length and width of the river is L x = 150 m , L y = 30 m L_x = 150 \, \text{m}, \, L_y = 30 \, \text{m} L x = 150 m , L y = 30 m .
The flow is assumed to be in the x x x -direction, i.e. u = 0.4 m s − 1 , v = 0 u = 0.4 \, \text{m} \, \text{s}^{−1}, v = 0 u = 0.4 m s − 1 , v = 0 .
The diffusion coefficients D x = D y = 0.5 m 2 s − 1 D_x = D_y = 0.5 \, \text{m}^2 \, \text{s}^{−1} D x = D y = 0.5 m 2 s − 1 .
The reaction rate K r = 0.01 s − 1 K_{\text{r}} = 0.01 \, \text{s}^{−1} K r = 0.01 s − 1 .
The mass transfer coefficient for oxygen K a = 0.02 s − 1 K_{\text{a}} = 0.02 \, \text{s}^{−1} K a = 0.02 s − 1 .
The saturation oxygen content O sat = 8 mg L − 1 O_{\text{sat}} = 8 \, \text{mg} \; \text{L}^{−1} O sat = 8 mg L − 1 .
The BOD discharge rate Q B = 70 mg s − 1 Q_{\text{B}} = 70 \, \text{mg} \; \text{s}^{−1} Q B = 70 mg s − 1 .
Place the point source Q B Q_{\text{B}} Q B at x = 10 m , y = 15 m x = 10 \, \text{m}, \, y = 15 \, \text{m} x = 10 m , y = 15 m .
Use Δ x = Δ y = 1 m \Delta x = \Delta y = 1 \, \text{m} Δ x = Δ y = 1 m for the homework. (You may use finer spacings).
5. Plot a x , y x, y x , y contour of the steady state (i.e. ∂ / ∂ t = 0 \partial / \partial t = 0 ∂ / ∂ t = 0 ) solution for both BOD and DO. (HINT: use Matlab tool: contour or contourf). Plot the DO sag curve and BOD curve along the center line y = 15 m y =15 \, \text{m} y = 15 m , i.e. plot B ( x , y = 15 m ) B(x, y = 15 \, \text{m}) B ( x , y = 15 m ) and O ( x , y = 15 m ) O(x, y = 15 \, \text{m}) O ( x , y = 15 m ) .
6. According to Environmental Protection Agency (EPA) requirement, the minimum steady state DO in a river must be above 6 mg L − 1 6 \, \text{mg} \; \text{L}^{−1} 6 mg L − 1 . Does the current discharge satisfy the EPA requirement? If not, propose alternative discharge methods, test with your model, and present your results. Note that you must maintain the total discharge rate at 70 mg s − 1 70 \, \text{mg} \; \text{s}^{−1} 70 mg s − 1 .
7. Suggest 3 or more ways to make your code more realisitic.
8. Suggest 3 or more ways to speed up your code.
本章讨论二维 Streeter-Phelps 模型 (1)(2) 的数值解. 用 forward in time, centered in space (FTCS) 方案, 将 Eq. (1) 离散为
1 τ ( B i , j n + 1 − B i , j n ) + u i , j n 2 h x ( B i + 1 , j n − B i − 1 , j n ) + v i , j n 2 h y ( B i , j + 1 n − B i , j − 1 n ) = D x h x 2 ( B i + 1 , j n − 2 B i , j n + B i − 1 , j n ) + D y h y 2 ( B i , j + 1 n − 2 B i , j n + B i , j − 1 n ) − K r B i , j n + Q B h x h y I i , j ( x 0 , y 0 ) , \begin{equation*}
\tag{3a}
\begin{aligned}
& \quad \; \frac{1}{\tau} \left( B^{n+1}_{i,j} - B^{n}_{i,j} \right) + \frac{u^{n}_{i, j}}{2 h_x} \left( B^{n}_{i+1,j} - B^{n}_{i-1,j} \right) + \frac{v^n_{i,j}}{2 h_y} \left( B^{n}_{i,j+1} - B^{n}_{i,j-1} \right) \\
& = \frac{D_x}{h^2_x} \left( B^{n}_{i+1,j} - 2 B^{n}_{i,j} + B^{n}_{i-1,j} \right) + \frac{D_y}{h^2_y} \left( B^{n}_{i,j+1} - 2 B^{n}_{i,j} + B^{n}_{i,j-1} \right) - K_{\text{r}} B^n_{i,j} + \frac{Q_\text{B}}{h_x h_y} \mathrm{I}_{i,j}(x_0, y_0),
\end{aligned}
\end{equation*} τ 1 ( B i , j n + 1 − B i , j n ) + 2 h x u i , j n ( B i + 1 , j n − B i − 1 , j n ) + 2 h y v i , j n ( B i , j + 1 n − B i , j − 1 n ) = h x 2 D x ( B i + 1 , j n − 2 B i , j n + B i − 1 , j n ) + h y 2 D y ( B i , j + 1 n − 2 B i , j n + B i , j − 1 n ) − K r B i , j n + h x h y Q B I i , j ( x 0 , y 0 ) , ( 3a )
where
n = 0 , 1 , ⋯ , N , i = 0 , 1 , ⋯ I , j = 0 , 1 , ⋯ , J , n = 0, 1, \cdots, N, \quad i = 0, 1, \cdots I, \quad j = 0, 1, \cdots, J, n = 0 , 1 , ⋯ , N , i = 0 , 1 , ⋯ I , j = 0 , 1 , ⋯ , J ,
示性函数
I i , j ( x 0 , y 0 ) : = { 1 , 源点 ( x 0 , y 0 ) 落在 ( i , j ) 格点所代表的那个有限体积元 , 0 , otherwise . \begin{equation*}
\tag{4}
\mathrm{I}_{i,j}(x_0, y_0) := \left\{
\begin{aligned}
& 1, \quad 源点 (x_0, y_0) \text{ 落在 $(i, j)$ 格点所代表的那个有限体积元}, \\
& 0, \quad \text{otherwise}.
\end{aligned}
\right.
\end{equation*} I i , j ( x 0 , y 0 ) := { 1 , 源点 ( x 0 , y 0 ) 落在 ( i , j ) 格点所代表的那个有限体积元 , 0 , otherwise . ( 4 )
可以设置多个源点. 将 Eq. (2) 离散为
1 τ ( O i , j n + 1 − O i , j n ) + u i , j n 2 h x ( O i + 1 , j n − O i − 1 , j n ) + v i , j n 2 h y ( O i , j + 1 n − O i , j − 1 n ) = D x h x 2 ( O i + 1 , j n − 2 O i , j n + O i − 1 , j n ) + D y h y 2 ( O i , j + 1 n − 2 O i , j n + O i , j − 1 n ) + K a ( O sat − O i , j n ) − K r B i , j n . \begin{equation*}
\tag{3b}
\begin{aligned}
& \quad \; \frac{1}{\tau} \left( O^{n+1}_{i,j} - O^{n}_{i,j} \right) + \frac{u^{n}_{i, j}}{2 h_x} \left( O^{n}_{i+1,j} - O^{n}_{i-1,j} \right) + \frac{v^n_{i,j}}{2 h_y} \left( O^{n}_{i,j+1} - O^{n}_{i,j-1} \right) \\
& = \frac{D_x}{h^2_x} \left( O^{n}_{i+1,j} - 2 O^{n}_{i,j} + O^{n}_{i-1,j} \right) + \frac{D_y}{h^2_y} \left( O^{n}_{i,j+1} - 2 O^{n}_{i,j} + O^{n}_{i,j-1} \right) + K_{\text{a}} \left( O_{\text{sat}} - O^n_{i,j} \right) - K_{\text{r}} B^n_{i,j}.
\end{aligned}
\end{equation*} τ 1 ( O i , j n + 1 − O i , j n ) + 2 h x u i , j n ( O i + 1 , j n − O i − 1 , j n ) + 2 h y v i , j n ( O i , j + 1 n − O i , j − 1 n ) = h x 2 D x ( O i + 1 , j n − 2 O i , j n + O i − 1 , j n ) + h y 2 D y ( O i , j + 1 n − 2 O i , j n + O i , j − 1 n ) + K a ( O sat − O i , j n ) − K r B i , j n . ( 3b )
格式 (3) 中的速度场 u n , v n u^n, v^n u n , v n 被假定已由数值模式的动力框架给出.
Problem 1.1
我们要在二维有界区域
t ≥ 0 , ( x , y ) ∈ Ω : = [ 0 , L x ] × [ 0 , L y ] \begin{equation*}
\tag{5}
t \ge 0, \quad (x, y) \in \Omega := [0, L_x] \times [0, L_y]
\end{equation*} t ≥ 0 , ( x , y ) ∈ Ω := [ 0 , L x ] × [ 0 , L y ] ( 5 )
上求解模型 (1)(2), 为此, 须为模型设置合适的初始条件和边界条件, 这里取为
{ ( I.C. ) F ∣ t = 0 = F 0 ( x , y ) , ( x , y ) ∈ Ω , ( B.C. ) { ( Dirichlet ) F ∣ x = 0 = F 0 ( y , t ) , ( Neumann ) F x ∣ x = L x = F 1 ( y , t ) , ( solid wall ) ( n ⋅ u ) y = 0 , L y = 0 , ( solid wall 2 ) [ n ⋅ ( − κ ∇ F ) ] y = 0 , L y = 0. \begin{equation*}
\tag{6}
\left\{
\begin{aligned}
& (\text{I.C.}) \quad
\begin{aligned}
% \tag{6}
\left. F \right|_{t = 0} = F^0(x, y), \quad (x,y) \in \Omega,
\end{aligned} \\
& (\text{B.C.}) \quad
\left\{
\begin{aligned}
% \tag{7}
& (\text{Dirichlet}) & \quad \left. F \right|_{x = 0} = F_0(y, t), \\
& (\text{Neumann}) & \quad \left. F_x \right|_{x = L_x} = F_1(y, t), \\
& (\text{solid wall}) & \quad \left( \boldsymbol{n} \cdot \boldsymbol{u} \right)_{y = 0, \, L_y} = 0, \\
& ({\color{blue}{\text{solid wall 2}}}) & \quad \left[ \boldsymbol{n} \cdot (- \kappa \nabla F) \right]_{y = 0, \, L_y} = 0.
\end{aligned}
\right. \\
\end{aligned}
\right.
\end{equation*} ⎩ ⎨ ⎧ ( I.C. ) F ∣ t = 0 = F 0 ( x , y ) , ( x , y ) ∈ Ω , ( B.C. ) ⎩ ⎨ ⎧ ( Dirichlet ) ( Neumann ) ( solid wall ) ( solid wall 2 ) F ∣ x = 0 = F 0 ( y , t ) , F x ∣ x = L x = F 1 ( y , t ) , ( n ⋅ u ) y = 0 , L y = 0 , [ n ⋅ ( − κ ∇ F ) ] y = 0 , L y = 0. ( 6 )
Eq. (6) 中, F F F 指代 B B B 和 O O O .
相应地, 为格式 (3) 引入由
F − 1 , j n + F 1 , j n 2 = F 0 ( t n , y j ) , F 1 , j n − F − 1 , j n 2 h x = F 1 ( t n , y j ) , F i , − 1 n = 3 F i , 0 n − 3 F i , 1 n + F i , 2 n , F i , J + 1 n = 3 F i , J n − 3 F i , J − 1 n + F i , J − 2 n , \begin{align*}
\tag{7a}
& \frac{F^n_{-1, j} + F^n_{1, j}}{2} = F_0(t_n, y_j), \\
\tag{7b}
& \frac{F^n_{1, j} - F^n_{-1, j}}{2 h_x} = F_1(t_n, y_j), \\
\tag{8a}
& \color{gray}{F^n_{i, -1} = 3 F^n_{i, 0} - 3F^n_{i, 1} + F^n_{i, 2}}, \\
\tag{8b}
& \color{gray}{F^n_{i, J+1} = 3 F^n_{i, J} - 3 F^n_{i, J-1} + F^n_{i, J-2}},
\end{align*} 2 F − 1 , j n + F 1 , j n = F 0 ( t n , y j ) , 2 h x F 1 , j n − F − 1 , j n = F 1 ( t n , y j ) , F i , − 1 n = 3 F i , 0 n − 3 F i , 1 n + F i , 2 n , F i , J + 1 n = 3 F i , J n − 3 F i , J − 1 n + F i , J − 2 n , ( 7a ) ( 7b ) ( 8a ) ( 8b )
决定的 halo points F − 1 , j n , F I + 1 , j n , F i , − 1 n , F i , J + 1 n F^n_{-1, j}, F^n_{I + 1, j}, F^n_{i, -1}, F^n_{i, J + 1} F − 1 , j n , F I + 1 , j n , F i , − 1 n , F i , J + 1 n . 其中, Eq. (8) 的设计考虑是使格式 (3) 在网格边界 j = 0 , J j = 0, J j = 0 , J 处能以 O ( h y ) \mathcal{O}(h_y) O ( h y ) 精度离散模型 (1)(2) 在边界 y = 0 , L y y = 0, L_y y = 0 , L y 处的 Fick 扩散项. Halo points F i , − 1 n , F i , J + 1 n F^n_{i, -1}, F^n_{i, J+1} F i , − 1 n , F i , J + 1 n 的设计 (8) 仅考虑 Fick 扩散, 是因为 (6) 中 B.C. 第二项 (y 边界 solid wall) 保证 v i , 0 n = v i , J n = 0 v^n_{i, 0} = v^n_{i, J} = 0 v i , 0 n = v i , J n = 0 , 所以无论 F i , − 1 n , F i , J + 1 n F^n_{i, -1}, F^n_{i, J+1} F i , − 1 n , F i , J + 1 n 如何选取, 都不会破坏格式 (3) 对模型 (1)(2) 在边界 y = 0 , L y y = 0, L_y y = 0 , L y 处 v v v -平流 (exactly zero) 的离散 . 一个重要问题是, 在物理上, solid wall 边界对 Fick 扩散的约束是怎样的? 若是梯度扩散通量为零, 即 n ⋅ ( − κ ∇ F ) = 0 \color{black}{ \boldsymbol{n} \cdot (- \kappa \nabla F) = 0} n ⋅ ( − κ ∇ F ) = 0 , 则与齐次 Neumann 边界 ∂ F / ∂ n = 0 \color{black}{\partial F / \partial n = 0} ∂ F / ∂ n = 0 等效? 发现 ChatGPT 支持 此观点, 且有符合常识的物理解释, 所以我将 (8) 改为
F i , − 1 n = F i , 1 n , F i , J + 1 n = F i , J − 1 n . \begin{equation*}
\tag{$8'$}
F^n_{i, -1} = F^n_{i, 1}, \quad F^n_{i, J+1} = F^n_{i, J-1}.
\end{equation*} F i , − 1 n = F i , 1 n , F i , J + 1 n = F i , J − 1 n . ( 8 ′ )
但, 若我坚持用 (8), 则会引入多大的误差呢? 毕竟, Eq. (8), 尽管不蕴含 (8'), 似乎也是在合理前提下 (∂ 2 F / ∂ n 2 \partial^2 F / \partial n^2 ∂ 2 F / ∂ n 2 连续到边界? 这个条件过于苛刻吗? 零通量边界会引起二阶导突变甚至不存在?) 导出的?
引入 halo points (7)(8') 后, 在速度场 u i , j n , v i , j n u^n_{i,j}, v^n_{i,j} u i , j n , v i , j n 已获得的前提下 (速度场的求解方案 is beyond the scope of this article), 只需根据模型 (1)(2) 的初始条件 (6) 制作初始场
B i , j 0 = B 0 ( x i , y j ) , O i , j 0 = O 0 ( x i , y j ) , x i = i h x , y j = j h y , \begin{equation*}
\tag{9}
B^0_{i,j} = B^0(x_i, y_j), \quad O^0_{i,j} = O^0(x_i, y_j), \quad x_i = i h_x, \, y_j = j h_y,
\end{equation*} B i , j 0 = B 0 ( x i , y j ) , O i , j 0 = O 0 ( x i , y j ) , x i = i h x , y j = j h y , ( 9 )
就可驱动格式 (3) 向前积分. 以后, 简单起见, 将速度场取为常向量, 即认为 u n ≡ u , v n ≡ v u^n \equiv u, \, v^n \equiv v u n ≡ u , v n ≡ v , where u , v u, v u , v 是已知的常数.
Problem 1.2
本节对格式 (3) 作 von Neumann 稳定性分析. 为简化问题, I ignored the BOD discharge term δ ( x − x 0 , y − y 0 ) Q B \delta(x −x_0, y − y_0) Q_{\text{B}} δ ( x − x 0 , y − y 0 ) Q B in Eq. (1), and the K a O sat − K r B K_{\text{a}} O_{\text{sat}} − K_{\text{r}} B K a O sat − K r B terms in Eq. (2). 简化后的格式 (3a) 和 (3b) 有相同形式的增长因子,
G = ( 1 − 2 r x − 2 r y − K τ ) + 2 r x cos ( ξ x h x ) − i C r x sin ( ξ x h x ) + 2 r y cos ( ξ y h y ) − i C r y sin ( ξ y h y ) , \begin{equation*}
\tag{10}
\begin{aligned}
G & = (1 - 2 r_x - 2 r_y - K \tau) \\
& + 2 r_x \cos{(\xi_x h_x)} - \mathrm{i} \, \mathrm{Cr}_x \, \sin{(\xi_x h_x)} \\
& + 2 r_y \cos{(\xi_y h_y)} - \mathrm{i} \, \mathrm{Cr}_y \, \sin{(\xi_y h_y)},
\end{aligned}
\end{equation*} G = ( 1 − 2 r x − 2 r y − K τ ) + 2 r x cos ( ξ x h x ) − i Cr x sin ( ξ x h x ) + 2 r y cos ( ξ y h y ) − i Cr y sin ( ξ y h y ) , ( 10 )
其中 K K K 指代 K r K_{\mathrm{r}} K r 和 K a K_{\mathrm{a}} K a ,
C r i : = u i τ / h i 2 , r i : = D τ / h i 2 , P e , i : = u i h / D , \begin{equation*}
\mathrm{Cr}_i := u_i \tau / h_i^2, \; r_i := D \tau / h_i^2, \; P_{\mathrm{e}, i} := u_i h / D,
\end{equation*} Cr i := u i τ / h i 2 , r i := D τ / h i 2 , P e , i := u i h / D ,
下标 i = 1 , 2 , 3 i = 1, 2, 3 i = 1 , 2 , 3 分别表示 x , y , z x, y, z x , y , z 方向. 要使格式稳定, 就是要选取 τ , h x , h y \tau, h_x, h_y τ , h x , h y , 使加强的 von Neumann 条件 (Hw 3 的 Eq. (1.22))
∣ G ∣ ≤ 1 , ∀ ∣ ξ x h x ∣ ≤ π , ∣ ξ y h y ∣ ≤ π \begin{equation*}
\tag{11}
|G| \le 1, \quad \forall |\xi_x h_x| \le \pi, \, |\xi_y h_y| \le \pi
\end{equation*} ∣ G ∣ ≤ 1 , ∀∣ ξ x h x ∣ ≤ π , ∣ ξ y h y ∣ ≤ π ( 11 )
被满足.
设 A A A 是定义了加法运算的集合, A 1 , A 2 ⊂ A A_1, A_2 \subset A A 1 , A 2 ⊂ A . 定义集合间的加法
A 1 + A 2 : = { a 1 + a 2 ∣ a i ∈ A i } . \begin{equation*}
\tag{12}
A_1 + A_2 := \left\{ \left. a_1 + a_2 \; \right| \; a_i \in A_i \right\}.
\end{equation*} A 1 + A 2 := { a 1 + a 2 ∣ a i ∈ A i } . ( 12 )
记平面直角坐标系上的椭圆
C ( a , b ) : = { ( x , y ) ∣ x 2 a 2 + y 2 b 2 = 1 } , a , b > 0. \begin{equation*}
\tag{13}
C(a, b) := \left\{(x, y) \; \left| \; \frac{x^2}{a^2} + \frac{y^2}{b^2} = 1 \right. \right\}, \quad a, b > 0.
\end{equation*} C ( a , b ) := { ( x , y ) a 2 x 2 + b 2 y 2 = 1 } , a , b > 0. ( 13 )
当 ( ξ x h x , ξ y h y ) (\xi_x h_x, \xi_y h_y) ( ξ x h x , ξ y h y ) 取遍 [ − π , π ) 2 [-\pi, \pi)^2 [ − π , π ) 2 , 由 (10) 知 G = x + i y G = x + \mathrm{i} y G = x + i y 在复平面上的轨迹是
z 0 + C ( 2 r x , ∣ C r x ∣ ) + C ( 2 r y , ∣ C r y ∣ ) , \begin{equation*}
\tag{14}
z_0 + C(2 r_x, |\mathrm{Cr}_x|) + C(2 r_y, |\mathrm{Cr}_y|),
\end{equation*} z 0 + C ( 2 r x , ∣ Cr x ∣ ) + C ( 2 r y , ∣ Cr y ∣ ) , ( 14 )
其中
z 0 : = 1 − 2 r x − 2 r y − K τ \begin{equation*}
\tag{15}
z_0 := 1 - 2 r_x - 2 r_y - K \tau
\end{equation*} z 0 := 1 − 2 r x − 2 r y − K τ ( 15 )
在实轴上. 由 (14) 看出, 命题 (11) 的一个充分条件是
∣ z 0 ∣ + max { 2 r x , ∣ C r x ∣ } + max { 2 r y , ∣ C r y ∣ } ≤ 1 , \begin{equation*}
\tag{16}
|z_0| + \max{\{2 r_x, |\mathrm{Cr}_x|\}} + \max{\{2 r_y, |\mathrm{Cr}_y|\}} \le 1,
\end{equation*} ∣ z 0 ∣ + max { 2 r x , ∣ Cr x ∣ } + max { 2 r y , ∣ Cr y ∣ } ≤ 1 , ( 16 )
一个必要条件是
∣ z 0 ∣ + 2 r x + 2 r y ≤ 1. \begin{equation*}
\tag{17}
|z_0| + 2 r_x + 2 r_y \le 1.
\end{equation*} ∣ z 0 ∣ + 2 r x + 2 r y ≤ 1. ( 17 )
故在
2 r x ≥ ∣ C r x ∣ , 2 r y ≥ ∣ C r y ∣ \begin{equation*}
\tag{18}
2 r_x \ge |\mathrm{Cr}_x|, \, 2 r_y \ge |\mathrm{Cr}_y|
\end{equation*} 2 r x ≥ ∣ Cr x ∣ , 2 r y ≥ ∣ Cr y ∣ ( 18 )
的大前提下, 命题 (11) 的充分条件 (16) 即是必要条件 (17), 从而是充要的. 对于 (18) 以外的情形, (16) 是充分但不必要的. 简单起见, 我们用命题 (16) 的、从而也是命题 (11) 的一个充分不必要 条件,
∣ Cr y ∣ ≤ ∣ Cr x ∣ ≤ 2 r x = 2 r y ≤ ( 1 − K τ ) / 2 , \begin{equation*}
\tag{19}
|\text{Cr}_y| \le |\text{Cr}_x| \le 2 r_x = 2 r_y \le (1 - K \tau) / 2,
\end{equation*} ∣ Cr y ∣ ≤ ∣ Cr x ∣ ≤ 2 r x = 2 r y ≤ ( 1 − K τ ) /2 , ( 19 )
来决定网格参数 h x , h y , τ h_x, h_y, \tau h x , h y , τ . 特别地, 当 v = 0 , D x = D y = D v = 0, \, D_x = D_y = D v = 0 , D x = D y = D , 条件 (19) 等价于
h x = h y = : h ≤ 2 D / ∣ u ∣ , τ ≤ ( 4 D / h 2 + max { K a , K r } ) − 1 . \begin{equation*}
h_x = h_y =: h \le {2 D / |u|}, \quad \tau \le \left(4 D / h^2 + \max{\{ K_{\mathrm{a}}, K_{\mathrm{r}} \}} \right)^{-1}.
\end{equation*} h x = h y =: h ≤ 2 D /∣ u ∣ , τ ≤ ( 4 D / h 2 + max { K a , K r } ) − 1 .
一般地, 要写出 (11) 的充要条件, 或许可先用 Lagrange multipliers (2019 , 2022 ) 方法写出 max ∣ ξ h ∣ ≤ π ∣ G ∣ \max_{|\xi h| \le \pi}|G| max ∣ ξ h ∣ ≤ π ∣ G ∣ 的表达式, 再由 (11) 获得 G G G 的其余参数应满足的约束. 计算似乎相当繁琐, 所以留给感兴趣的读者. 我用 Mathematica (MMA ) 软件算得 (11) 的一个充要条件是
MMA 似乎算不出来? 是我给的命令不好吗? \begin{equation*}
\tag{20}
\text{MMA 似乎算不出来? 是我给的命令不好吗?}
\end{equation*} MMA 似乎算不出来 ? 是我给的命令不好吗 ? ( 20 )
G = 1 - 2 rx - 2 ry - k tau + 2 rx cx - I crx sx + 2 ry cy - I crx sy
ForAll[ { cx, sx, cy, sy} , cx^ 2 + sx^ 2 == 1 & & cy^ 2 + sy^ 2 == 1 , Abs [ G] <= 1 ]
若作简化
Cr y = 0 , r x = r y = r , K τ = 0 , \begin{equation*}
\tag{21}
\text{Cr}_y = 0, \; r_x = r_y = r, \; K\tau = 0,
\end{equation*} Cr y = 0 , r x = r y = r , K τ = 0 , ( 21 )
也就是将命题 (21) 添加到大前提, 则用 MMA 显式化命题 (11)
Resolve[ ForAll[ { cx, sx, cy} , cx^ 2 + sx^ 2 == 1 & & - 1 <= cy <= 1 , ( crx* sx) ^ 2 + ( 2 * cx* r + 2 * cy* r - 2 * r - 2 * r + 1 ) ^ 2 <= 1 ] , Reals]
FullSimplify[ %, r >= 0 ]
将得到形式很复杂的复合命题
crx 2 + 4 ( r − 1 ) r ≤ 0 ∧ 4 r ≤ 1 ∧ ( crx = 2 r ∨ r = 0 ∨ crx + 2 r = 0 ∨ ( ( ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≤ 0 ∨ ( r 2 ( 2 r − 1 ) ( 4 r 2 − crx 2 ) ≤ 0 ∧ ( crx 2 + 4 r 2 ≠ 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) ≠ 0 ) ) ∨ 6 r < 1 ∨ 2 r 2 ≥ r ∨ r 2 ≤ 0 ) ∧ ( crx ≠ 0 ∨ crx 2 r ( 4 r − 1 ) ≥ 0 ∨ crx 2 r ( 4 r − 1 ) ( crx 4 − 4 crx 2 r ( r + 1 ) + 4 r 2 ) ≥ 0 ∨ ( r 2 ( crx 10 − 4 crx 8 ( r ( 9 r − 4 ) + 1 ) + 8 crx 6 r 2 ( 62 r 2 − 48 r + 11 ) + 16 crx 4 r 2 ( 4 r ( r ( 3 r ( 3 r − 4 ) + 10 ) − 3 ) + 1 ) − 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) + 64 r 6 ) ≥ 0 ∧ ( crx 2 + 4 r 2 ≠ 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) ≠ 0 ) ) ∨ ( r 2 ( crx 10 + 4 crx 8 r ( 7 r − 4 ) − 8 crx 6 r 2 ( 2 r ( r + 8 ) − 7 ) + 16 crx 4 r 2 ( 4 r ( r ( 5 r ( 5 r − 4 ) + 11 ) − 3 ) + 1 ) − 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) + 64 r 6 ) ≤ 0 ∧ ( crx 2 + 4 r 2 ≠ 0 ∨ crx 4 + 4 crx 2 r ( 3 r − 2 ) + 8 r 2 ( 8 r 2 − 4 r + 1 ) ≠ 0 ) ∧ ( crx 10 + 4 crx 8 r ( 7 r − 4 ) + 16 crx 4 r 2 ( 4 r ( r ( 5 r ( 5 r − 4 ) + 11 ) − 3 ) + 1 ) + 64 r 6 ≠ 8 crx 6 r 2 ( 2 r ( r + 8 ) − 7 ) + 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) ∨ crx 2 r 2 + 4 r 4 ≤ 0 ∨ ( crx 2 + 4 r 2 ≠ 0 ∧ crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 ≠ 0 ) ) ) ∨ ( crx 4 + 4 crx 2 r ( 7 r − 3 ) + 4 r 2 ≠ 0 ∧ crx 2 r ( 4 r − 1 ) ( crx 4 + 4 crx 2 r ( 7 r − 3 ) + 4 r 2 ) ≤ 0 ) ∨ ( crx 2 r ( 4 r − 1 ) ( crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 ) ≤ 0 ∧ ( r 2 ≤ 0 ∨ crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 ≠ 0 ) ) ) ∧ ( crx = 0 ∨ ( crx 2 = 2 r ∧ r > 0 ) ∨ ( ( crx 2 + 4 r 2 ≠ 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) ≠ 0 ) ∧ crx 8 r 2 + 16 crx 2 r 6 ≤ 4 crx 4 ( crx 2 + 1 ) r 4 ) ∨ r = 0 ∨ ( crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ≠ 0 ∧ crx 2 r 2 ( crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ) ≥ 0 ∧ crx 2 r ( 6 r − 1 ) ≤ 0 ) ∨ crx 2 r ( 2 r − 1 ) ≥ 0 ∨ crx 2 ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≥ 0 ) ∧ ( 2 r 2 ≥ r ∨ ( crx 2 − 2 r ) 2 ≤ 0 ∨ r ( 2 r − 1 ) ( crx 2 − 2 r ) > 0 ) ∧ ( ( r > 0 ∧ ( r ( 2 r − 1 ) ( crx 2 − 2 r ) < 0 ∨ ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) > 0 ∧ 6 r < 1 ) ) ) ∨ crx 2 + 4 r 2 = 0 ∨ r 2 ( 2 r − 1 ) ( 16 r 4 − crx 4 ) ≤ 0 ∨ crx 8 r 2 + 64 r 8 ≤ 4 crx 4 ( 4 r 6 + r 4 ) ) ) ) ∧ ( crx = 2 r ∨ r = 0 ∨ crx + 2 r = 0 ∨ ( ( crx ≠ 0 ∨ crx 2 r ( 4 r − 1 ) ≥ 0 ∨ crx 2 r ( 4 r − 1 ) ( crx 4 − 4 crx 2 r ( r + 1 ) + 4 r 2 ) ≥ 0 ∨ ( r 2 ( crx 10 − 4 crx 8 ( r ( 9 r − 4 ) + 1 ) + 8 crx 6 r 2 ( 62 r 2 − 48 r + 11 ) + 16 crx 4 r 2 ( 4 r ( r ( 3 r ( 3 r − 4 ) + 10 ) − 3 ) + 1 ) − 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) + 64 r 6 ) ≥ 0 ∧ ( crx 2 + 4 r 2 ≠ 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) ≠ 0 ) ) ∨ ( r 2 ( crx 10 + 4 crx 8 r ( 7 r − 4 ) − 8 crx 6 r 2 ( 2 r ( r + 8 ) − 7 ) + 16 crx 4 r 2 ( 4 r ( r ( 5 r ( 5 r − 4 ) + 11 ) − 3 ) + 1 ) − 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) + 64 r 6 ) ≤ 0 ∧ ( crx 2 + 4 r 2 ≠ 0 ∨ crx 4 + 4 crx 2 r ( 3 r − 2 ) + 8 r 2 ( 8 r 2 − 4 r + 1 ) ≠ 0 ) ∧ ( crx 10 + 4 crx 8 r ( 7 r − 4 ) + 16 crx 4 r 2 ( 4 r ( r ( 5 r ( 5 r − 4 ) + 11 ) − 3 ) + 1 ) + 64 r 6 ≠ 8 crx 6 r 2 ( 2 r ( r + 8 ) − 7 ) + 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) ∨ crx 2 r 2 + 4 r 4 ≤ 0 ∨ ( crx 2 + 4 r 2 ≠ 0 ∧ crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 ≠ 0 ) ) ) ∨ ( crx 4 + 4 crx 2 r ( 7 r − 3 ) + 4 r 2 ≠ 0 ∧ crx 2 r ( 4 r − 1 ) ( crx 4 + 4 crx 2 r ( 7 r − 3 ) + 4 r 2 ) ≤ 0 ) ∨ crx 2 r ( 4 r − 1 ) ( crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 ) ≥ 0 ) ∧ ( crx = 0 ∨ ( crx 2 = 2 r ∧ r > 0 ) ∨ ( ( crx 2 + 4 r 2 ≠ 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) ≠ 0 ) ∧ crx 8 r 2 + 16 crx 2 r 6 ≤ 4 crx 4 ( crx 2 + 1 ) r 4 ) ∨ ( crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ≠ 0 ∧ crx 2 r 2 ( crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ) ≥ 0 ∧ crx 2 r ( 6 r − 1 ) ≤ 0 ) ∨ crx 2 r ( 2 r − 1 ) ≥ 0 ∨ crx 2 ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≥ 0 ∨ crx 2 r 2 > 0 ) ∧ ( 6 r 2 ≤ r ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ≤ 0 ∨ r ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) ) > 0 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 8 r 2 ( 6 r ( 3 r − 1 ) + 1 ) < 0 ) ∧ ( ( r > 0 ∧ ( r ( 2 r − 1 ) ( crx 2 − 2 r ) < 0 ∨ ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) > 0 ∧ 6 r < 1 ) ) ) ∨ crx 2 + 4 r 2 = 0 ∨ r 2 ( 2 r − 1 ) ( 16 r 4 − crx 4 ) ≤ 0 ∨ crx 2 r 2 + 4 r 4 > 0 ∨ crx 8 r 2 + 64 r 8 ≤ 4 crx 4 ( 4 r 6 + r 4 ) ) ) ) ∧ ( r = 0 ∨ crx 2 ≥ 4 r 2 ∨ ( ( crx = 0 ∨ crx 2 ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≥ 0 ∨ ( ( crx 2 r ( 6 r − 1 ) ≤ 0 ∨ ( r > 0 ∧ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ≤ 0 ) ) ∧ ( crx 2 r ( 6 r − 1 ) < 0 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ≠ 0 ) ) ∨ crx 2 r ( 2 r − 1 ) ≥ 0 ∨ ( crx 8 r 2 + 16 crx 2 r 6 ≤ 4 crx 4 ( crx 2 + 1 ) r 4 ∧ ( crx 2 + 4 r 2 ≠ 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) ≠ 0 ) ∧ ( crx 2 ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ( crx 2 + 4 r 2 ) 2 ≥ 0 ∨ ( crx ≠ 2 r ∧ crx 2 ≠ 2 r ∧ crx + 2 r ≠ 0 ∧ crx 2 + 2 r ≠ 0 ) ) ) ) ∧ ( 6 r = 1 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ≥ 0 ∨ r ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) ) > 0 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 8 r 2 ( 6 r ( 3 r − 1 ) + 1 ) < 0 ) ∧ ( crx = 2 r ∨ ( 0 < r < 1 6 ∧ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 8 r 2 ( 6 r ( 3 r − 1 ) + 1 ) > 0 ) ∨ crx + 2 r = 0 ∨ crx 2 + 4 r 2 = 0 ∨ ( crx 2 − 1 ) r 2 ( crx 4 − 16 r 4 ) ≥ 0 ) ∧ ( crx = 2 r ∨ 0 < r < 1 6 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 8 r 2 ( 6 r ( 3 r − 1 ) + 1 ) < 0 ∨ crx + 2 r = 0 ∨ crx 2 + 4 r 2 = 0 ∨ ( crx 2 − 1 ) r 2 ( crx 4 − 16 r 4 ) ≥ 0 ) ∧ ( crx = 2 r ∨ r ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) ) < 0 ∨ crx + 2 r = 0 ∨ crx 2 + 4 r 2 = 0 ∨ 0 < r < 1 6 ∨ r 2 ( 2 r − 1 ) ( 16 r 4 − crx 4 ) < 0 ∨ 4 crx 4 ( 4 r 6 + r 4 ) ≤ crx 8 r 2 + 64 r 8 ) ∧ ( crx = 2 r ∨ r ( 2 r − 1 ) ( crx 2 − 2 r ) < 0 ∨ crx + 2 r = 0 ∨ crx 2 + 4 r 2 = 0 ∨ ( crx 2 ( 3 − 10 r ) < 2 r ( 2 r ( 8 r − 3 ) + 1 ) ∧ 6 r < 1 ) ∨ r 2 ( 2 r − 1 ) ( 16 r 4 − crx 4 ) < 0 ∨ 4 crx 4 ( 4 r 6 + r 4 ) ≤ crx 8 r 2 + 64 r 8 ) ∧ ( 6 r < 1 ∨ 2 r 2 ≥ r ∨ ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≥ 0 ∨ ( ( crx 2 + 4 r 2 ≠ 0 ∨ crx 4 + 4 crx 2 r ( 3 r − 2 ) + 8 r 2 ( 8 r 2 − 4 r + 1 ) ≠ 0 ) ∧ ( ( crx ≠ 2 r ∧ crx + 2 r ≠ 0 ) ∨ crx 2 r 2 + 4 r 4 ≤ 0 ) ) ∨ ( r 2 ( 2 r − 1 ) ( 4 r 2 − crx 2 ) ≤ 0 ∧ ( crx 2 + 4 r 2 ≠ 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) ≠ 0 ) ) ) ) ) \begin{equation*}
\tag{22}
\text{crx}^2+4 (r-1) r\leq 0\land 4 r\leq 1\land \left(\text{crx}=2 r\lor r=0\lor \text{crx}+2 r=0\lor \left(\left(\left(\text{crx}^2-1\right) r^2 \left(\text{crx}^2-4 r^2\right)\leq 0\lor \left(r^2 (2 r-1) \left(4 r^2-\text{crx}^2\right)\leq 0\land \left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^2 (5 r-2)+2 r \left(6 r^2-4 r+1\right)\neq 0\right)\right)\lor 6 r<1\lor 2 r^2\geq r\lor r^2\leq 0\right)\land \left(\text{crx}\neq 0\lor \text{crx}^2 r (4 r-1)\geq 0\lor \text{crx}^2 r (4 r-1) \left(\text{crx}^4-4 \text{crx}^2 r (r+1)+4 r^2\right)\geq 0\lor \left(r^2 \left(\text{crx}^{10}-4 \text{crx}^8 (r (9 r-4)+1)+8 \text{crx}^6 r^2 \left(62 r^2-48 r+11\right)+16 \text{crx}^4 r^2 (4 r (r (3 r (3 r-4)+10)-3)+1)-16 \text{crx}^2 r^4 (8 r (5 r-2)+3)+64 r^6\right)\geq 0\land \left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^2 (5 r-2)+2 r \left(6 r^2-4 r+1\right)\neq 0\right)\right)\lor \left(r^2 \left(\text{crx}^{10}+4 \text{crx}^8 r (7 r-4)-8 \text{crx}^6 r^2 (2 r (r+8)-7)+16 \text{crx}^4 r^2 (4 r (r (5 r (5 r-4)+11)-3)+1)-16 \text{crx}^2 r^4 (8 r (5 r-2)+3)+64 r^6\right)\leq 0\land \left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^4+4 \text{crx}^2 r (3 r-2)+8 r^2 \left(8 r^2-4 r+1\right)\neq 0\right)\land \left(\text{crx}^{10}+4 \text{crx}^8 r (7 r-4)+16 \text{crx}^4 r^2 (4 r (r (5 r (5 r-4)+11)-3)+1)+64 r^6\neq 8 \text{crx}^6 r^2 (2 r (r+8)-7)+16 \text{crx}^2 r^4 (8 r (5 r-2)+3)\lor \text{crx}^2 r^2+4 r^4\leq 0\lor \left(\text{crx}^2+4 r^2\neq 0\land \text{crx}^4+\text{crx}^2 (4 (2-5 r) r-2)+4 r^2\neq 0\right)\right)\right)\lor \left(\text{crx}^4+4 \text{crx}^2 r (7 r-3)+4 r^2\neq 0\land \text{crx}^2 r (4 r-1) \left(\text{crx}^4+4 \text{crx}^2 r (7 r-3)+4 r^2\right)\leq 0\right)\lor \left(\text{crx}^2 r (4 r-1) \left(\text{crx}^4+\text{crx}^2 (4 (2-5 r) r-2)+4 r^2\right)\leq 0\land \left(r^2\leq 0\lor \text{crx}^4+\text{crx}^2 (4 (2-5 r) r-2)+4 r^2\neq 0\right)\right)\right)\land \left(\text{crx}=0\lor \left(\text{crx}^2=2 r\land r>0\right)\lor \left(\left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^2 (5 r-2)+2 r \left(6 r^2-4 r+1\right)\neq 0\right)\land \text{crx}^8 r^2+16 \text{crx}^2 r^6\leq 4 \text{crx}^4 \left(\text{crx}^2+1\right) r^4\right)\lor r=0\lor \left(\text{crx}^4+4 \text{crx}^2 r (8 r-3)+4 r^2\neq 0\land \text{crx}^2 r^2 \left(\text{crx}^4+4 \text{crx}^2 r (8 r-3)+4 r^2\right)\geq 0\land \text{crx}^2 r (6 r-1)\leq 0\right)\lor \text{crx}^2 r (2 r-1)\geq 0\lor \text{crx}^2 \left(\text{crx}^2-1\right) r^2 \left(\text{crx}^2-4 r^2\right)\geq 0\right)\land \left(2 r^2\geq r\lor \left(\text{crx}^2-2 r\right)^2\leq 0\lor r (2 r-1) \left(\text{crx}^2-2 r\right)>0\right)\land \left(\left(r>0\land \left(r (2 r-1) \left(\text{crx}^2-2 r\right)<0\lor \left(\text{crx}^2 (10 r-3)+2 r (2 r (8 r-3)+1)>0\land 6 r<1\right)\right)\right)\lor \text{crx}^2+4 r^2=0\lor r^2 (2 r-1) \left(16 r^4-\text{crx}^4\right)\leq 0\lor \text{crx}^8 r^2+64 r^8\leq 4 \text{crx}^4 \left(4 r^6+r^4\right)\right)\right)\right)\land \left(\text{crx}=2 r\lor r=0\lor \text{crx}+2 r=0\lor \left(\left(\text{crx}\neq 0\lor \text{crx}^2 r (4 r-1)\geq 0\lor \text{crx}^2 r (4 r-1) \left(\text{crx}^4-4 \text{crx}^2 r (r+1)+4 r^2\right)\geq 0\lor \left(r^2 \left(\text{crx}^{10}-4 \text{crx}^8 (r (9 r-4)+1)+8 \text{crx}^6 r^2 \left(62 r^2-48 r+11\right)+16 \text{crx}^4 r^2 (4 r (r (3 r (3 r-4)+10)-3)+1)-16 \text{crx}^2 r^4 (8 r (5 r-2)+3)+64 r^6\right)\geq 0\land \left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^2 (5 r-2)+2 r \left(6 r^2-4 r+1\right)\neq 0\right)\right)\lor \left(r^2 \left(\text{crx}^{10}+4 \text{crx}^8 r (7 r-4)-8 \text{crx}^6 r^2 (2 r (r+8)-7)+16 \text{crx}^4 r^2 (4 r (r (5 r (5 r-4)+11)-3)+1)-16 \text{crx}^2 r^4 (8 r (5 r-2)+3)+64 r^6\right)\leq 0\land \left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^4+4 \text{crx}^2 r (3 r-2)+8 r^2 \left(8 r^2-4 r+1\right)\neq 0\right)\land \left(\text{crx}^{10}+4 \text{crx}^8 r (7 r-4)+16 \text{crx}^4 r^2 (4 r (r (5 r (5 r-4)+11)-3)+1)+64 r^6\neq 8 \text{crx}^6 r^2 (2 r (r+8)-7)+16 \text{crx}^2 r^4 (8 r (5 r-2)+3)\lor \text{crx}^2 r^2+4 r^4\leq 0\lor \left(\text{crx}^2+4 r^2\neq 0\land \text{crx}^4+\text{crx}^2 (4 (2-5 r) r-2)+4 r^2\neq 0\right)\right)\right)\lor \left(\text{crx}^4+4 \text{crx}^2 r (7 r-3)+4 r^2\neq 0\land \text{crx}^2 r (4 r-1) \left(\text{crx}^4+4 \text{crx}^2 r (7 r-3)+4 r^2\right)\leq 0\right)\lor \text{crx}^2 r (4 r-1) \left(\text{crx}^4+\text{crx}^2 (4 (2-5 r) r-2)+4 r^2\right)\geq 0\right)\land \left(\text{crx}=0\lor \left(\text{crx}^2=2 r\land r>0\right)\lor \left(\left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^2 (5 r-2)+2 r \left(6 r^2-4 r+1\right)\neq 0\right)\land \text{crx}^8 r^2+16 \text{crx}^2 r^6\leq 4 \text{crx}^4 \left(\text{crx}^2+1\right) r^4\right)\lor \left(\text{crx}^4+4 \text{crx}^2 r (8 r-3)+4 r^2\neq 0\land \text{crx}^2 r^2 \left(\text{crx}^4+4 \text{crx}^2 r (8 r-3)+4 r^2\right)\geq 0\land \text{crx}^2 r (6 r-1)\leq 0\right)\lor \text{crx}^2 r (2 r-1)\geq 0\lor \text{crx}^2 \left(\text{crx}^2-1\right) r^2 \left(\text{crx}^2-4 r^2\right)\geq 0\lor \text{crx}^2 r^2>0\right)\land \left(6 r^2\leq r\lor \text{crx}^4+4 \text{crx}^2 r (8 r-3)+4 r^2\leq 0\lor r \left(\text{crx}^2 (10 r-3)+2 r (2 r (8 r-3)+1)\right)>0\lor \text{crx}^4+4 \text{crx}^2 r (8 r-3)+8 r^2 (6 r (3 r-1)+1)<0\right)\land \left(\left(r>0\land \left(r (2 r-1) \left(\text{crx}^2-2 r\right)<0\lor \left(\text{crx}^2 (10 r-3)+2 r (2 r (8 r-3)+1)>0\land 6 r<1\right)\right)\right)\lor \text{crx}^2+4 r^2=0\lor r^2 (2 r-1) \left(16 r^4-\text{crx}^4\right)\leq 0\lor \text{crx}^2 r^2+4 r^4>0\lor \text{crx}^8 r^2+64 r^8\leq 4 \text{crx}^4 \left(4 r^6+r^4\right)\right)\right)\right)\land \left(r=0\lor \text{crx}^2\geq 4 r^2\lor \left(\left(\text{crx}=0\lor \text{crx}^2 \left(\text{crx}^2-1\right) r^2 \left(\text{crx}^2-4 r^2\right)\geq 0\lor \left(\left(\text{crx}^2 r (6 r-1)\leq 0\lor \left(r>0\land \text{crx}^4+4 \text{crx}^2 r (8 r-3)+4 r^2\leq 0\right)\right)\land \left(\text{crx}^2 r (6 r-1)<0\lor \text{crx}^4+4 \text{crx}^2 r (8 r-3)+4 r^2\neq 0\right)\right)\lor \text{crx}^2 r (2 r-1)\geq 0\lor \left(\text{crx}^8 r^2+16 \text{crx}^2 r^6\leq 4 \text{crx}^4 \left(\text{crx}^2+1\right) r^4\land \left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^2 (5 r-2)+2 r \left(6 r^2-4 r+1\right)\neq 0\right)\land \left(\text{crx}^2 \left(\text{crx}^2-1\right) r^2 \left(\text{crx}^2-4 r^2\right) \left(\text{crx}^2+4 r^2\right)^2\geq 0\lor \left(\text{crx}\neq 2 r\land \text{crx}^2\neq 2 r\land \text{crx}+2 r\neq 0\land \text{crx}^2+2 r\neq 0\right)\right)\right)\right)\land \left(6 r=1\lor \text{crx}^4+4 \text{crx}^2 r (8 r-3)+4 r^2\geq 0\lor r \left(\text{crx}^2 (10 r-3)+2 r (2 r (8 r-3)+1)\right)>0\lor \text{crx}^4+4 \text{crx}^2 r (8 r-3)+8 r^2 (6 r (3 r-1)+1)<0\right)\land \left(\text{crx}=2 r\lor \left(0<r<\frac{1}{6}\land \text{crx}^4+4 \text{crx}^2 r (8 r-3)+8 r^2 (6 r (3 r-1)+1)>0\right)\lor \text{crx}+2 r=0\lor \text{crx}^2+4 r^2=0\lor \left(\text{crx}^2-1\right) r^2 \left(\text{crx}^4-16 r^4\right)\geq 0\right)\land \left(\text{crx}=2 r\lor 0<r<\frac{1}{6}\lor \text{crx}^4+4 \text{crx}^2 r (8 r-3)+8 r^2 (6 r (3 r-1)+1)<0\lor \text{crx}+2 r=0\lor \text{crx}^2+4 r^2=0\lor \left(\text{crx}^2-1\right) r^2 \left(\text{crx}^4-16 r^4\right)\geq 0\right)\land \left(\text{crx}=2 r\lor r \left(\text{crx}^2 (10 r-3)+2 r (2 r (8 r-3)+1)\right)<0\lor \text{crx}+2 r=0\lor \text{crx}^2+4 r^2=0\lor 0<r<\frac{1}{6}\lor r^2 (2 r-1) \left(16 r^4-\text{crx}^4\right)<0\lor 4 \text{crx}^4 \left(4 r^6+r^4\right)\leq \text{crx}^8 r^2+64 r^8\right)\land \left(\text{crx}=2 r\lor r (2 r-1) \left(\text{crx}^2-2 r\right)<0\lor \text{crx}+2 r=0\lor \text{crx}^2+4 r^2=0\lor \left(\text{crx}^2 (3-10 r)<2 r (2 r (8 r-3)+1)\land 6 r<1\right)\lor r^2 (2 r-1) \left(16 r^4-\text{crx}^4\right)<0\lor 4 \text{crx}^4 \left(4 r^6+r^4\right)\leq \text{crx}^8 r^2+64 r^8\right)\land \left(6 r<1\lor 2 r^2\geq r\lor \left(\text{crx}^2-1\right) r^2 \left(\text{crx}^2-4 r^2\right)\geq 0\lor \left(\left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^4+4 \text{crx}^2 r (3 r-2)+8 r^2 \left(8 r^2-4 r+1\right)\neq 0\right)\land \left((\text{crx}\neq 2 r\land \text{crx}+2 r\neq 0)\lor \text{crx}^2 r^2+4 r^4\leq 0\right)\right)\lor \left(r^2 (2 r-1) \left(4 r^2-\text{crx}^2\right)\leq 0\land \left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^2 (5 r-2)+2 r \left(6 r^2-4 r+1\right)\neq 0\right)\right)\right)\right)\right)
\end{equation*} crx 2 + 4 ( r − 1 ) r ≤ 0 ∧ 4 r ≤ 1 ∧ ( crx = 2 r ∨ r = 0 ∨ crx + 2 r = 0 ∨ ( ( ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≤ 0 ∨ ( r 2 ( 2 r − 1 ) ( 4 r 2 − crx 2 ) ≤ 0 ∧ ( crx 2 + 4 r 2 = 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) = 0 ) ) ∨ 6 r < 1 ∨ 2 r 2 ≥ r ∨ r 2 ≤ 0 ) ∧ ( crx = 0 ∨ crx 2 r ( 4 r − 1 ) ≥ 0 ∨ crx 2 r ( 4 r − 1 ) ( crx 4 − 4 crx 2 r ( r + 1 ) + 4 r 2 ) ≥ 0 ∨ ( r 2 ( crx 10 − 4 crx 8 ( r ( 9 r − 4 ) + 1 ) + 8 crx 6 r 2 ( 62 r 2 − 48 r + 11 ) + 16 crx 4 r 2 ( 4 r ( r ( 3 r ( 3 r − 4 ) + 10 ) − 3 ) + 1 ) − 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) + 64 r 6 ) ≥ 0 ∧ ( crx 2 + 4 r 2 = 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) = 0 ) ) ∨ ( r 2 ( crx 10 + 4 crx 8 r ( 7 r − 4 ) − 8 crx 6 r 2 ( 2 r ( r + 8 ) − 7 ) + 16 crx 4 r 2 ( 4 r ( r ( 5 r ( 5 r − 4 ) + 11 ) − 3 ) + 1 ) − 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) + 64 r 6 ) ≤ 0 ∧ ( crx 2 + 4 r 2 = 0 ∨ crx 4 + 4 crx 2 r ( 3 r − 2 ) + 8 r 2 ( 8 r 2 − 4 r + 1 ) = 0 ) ∧ ( crx 10 + 4 crx 8 r ( 7 r − 4 ) + 16 crx 4 r 2 ( 4 r ( r ( 5 r ( 5 r − 4 ) + 11 ) − 3 ) + 1 ) + 64 r 6 = 8 crx 6 r 2 ( 2 r ( r + 8 ) − 7 ) + 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) ∨ crx 2 r 2 + 4 r 4 ≤ 0 ∨ ( crx 2 + 4 r 2 = 0 ∧ crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 = 0 ) ) ) ∨ ( crx 4 + 4 crx 2 r ( 7 r − 3 ) + 4 r 2 = 0 ∧ crx 2 r ( 4 r − 1 ) ( crx 4 + 4 crx 2 r ( 7 r − 3 ) + 4 r 2 ) ≤ 0 ) ∨ ( crx 2 r ( 4 r − 1 ) ( crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 ) ≤ 0 ∧ ( r 2 ≤ 0 ∨ crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 = 0 ) ) ) ∧ ( crx = 0 ∨ ( crx 2 = 2 r ∧ r > 0 ) ∨ ( ( crx 2 + 4 r 2 = 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) = 0 ) ∧ crx 8 r 2 + 16 crx 2 r 6 ≤ 4 crx 4 ( crx 2 + 1 ) r 4 ) ∨ r = 0 ∨ ( crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 = 0 ∧ crx 2 r 2 ( crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ) ≥ 0 ∧ crx 2 r ( 6 r − 1 ) ≤ 0 ) ∨ crx 2 r ( 2 r − 1 ) ≥ 0 ∨ crx 2 ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≥ 0 ) ∧ ( 2 r 2 ≥ r ∨ ( crx 2 − 2 r ) 2 ≤ 0 ∨ r ( 2 r − 1 ) ( crx 2 − 2 r ) > 0 ) ∧ ( ( r > 0 ∧ ( r ( 2 r − 1 ) ( crx 2 − 2 r ) < 0 ∨ ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) > 0 ∧ 6 r < 1 ) ) ) ∨ crx 2 + 4 r 2 = 0 ∨ r 2 ( 2 r − 1 ) ( 16 r 4 − crx 4 ) ≤ 0 ∨ crx 8 r 2 + 64 r 8 ≤ 4 crx 4 ( 4 r 6 + r 4 ) ) ) ) ∧ ( crx = 2 r ∨ r = 0 ∨ crx + 2 r = 0 ∨ ( ( crx = 0 ∨ crx 2 r ( 4 r − 1 ) ≥ 0 ∨ crx 2 r ( 4 r − 1 ) ( crx 4 − 4 crx 2 r ( r + 1 ) + 4 r 2 ) ≥ 0 ∨ ( r 2 ( crx 10 − 4 crx 8 ( r ( 9 r − 4 ) + 1 ) + 8 crx 6 r 2 ( 62 r 2 − 48 r + 11 ) + 16 crx 4 r 2 ( 4 r ( r ( 3 r ( 3 r − 4 ) + 10 ) − 3 ) + 1 ) − 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) + 64 r 6 ) ≥ 0 ∧ ( crx 2 + 4 r 2 = 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) = 0 ) ) ∨ ( r 2 ( crx 10 + 4 crx 8 r ( 7 r − 4 ) − 8 crx 6 r 2 ( 2 r ( r + 8 ) − 7 ) + 16 crx 4 r 2 ( 4 r ( r ( 5 r ( 5 r − 4 ) + 11 ) − 3 ) + 1 ) − 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) + 64 r 6 ) ≤ 0 ∧ ( crx 2 + 4 r 2 = 0 ∨ crx 4 + 4 crx 2 r ( 3 r − 2 ) + 8 r 2 ( 8 r 2 − 4 r + 1 ) = 0 ) ∧ ( crx 10 + 4 crx 8 r ( 7 r − 4 ) + 16 crx 4 r 2 ( 4 r ( r ( 5 r ( 5 r − 4 ) + 11 ) − 3 ) + 1 ) + 64 r 6 = 8 crx 6 r 2 ( 2 r ( r + 8 ) − 7 ) + 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) ∨ crx 2 r 2 + 4 r 4 ≤ 0 ∨ ( crx 2 + 4 r 2 = 0 ∧ crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 = 0 ) ) ) ∨ ( crx 4 + 4 crx 2 r ( 7 r − 3 ) + 4 r 2 = 0 ∧ crx 2 r ( 4 r − 1 ) ( crx 4 + 4 crx 2 r ( 7 r − 3 ) + 4 r 2 ) ≤ 0 ) ∨ crx 2 r ( 4 r − 1 ) ( crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 ) ≥ 0 ) ∧ ( crx = 0 ∨ ( crx 2 = 2 r ∧ r > 0 ) ∨ ( ( crx 2 + 4 r 2 = 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) = 0 ) ∧ crx 8 r 2 + 16 crx 2 r 6 ≤ 4 crx 4 ( crx 2 + 1 ) r 4 ) ∨ ( crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 = 0 ∧ crx 2 r 2 ( crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ) ≥ 0 ∧ crx 2 r ( 6 r − 1 ) ≤ 0 ) ∨ crx 2 r ( 2 r − 1 ) ≥ 0 ∨ crx 2 ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≥ 0 ∨ crx 2 r 2 > 0 ) ∧ ( 6 r 2 ≤ r ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ≤ 0 ∨ r ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) ) > 0 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 8 r 2 ( 6 r ( 3 r − 1 ) + 1 ) < 0 ) ∧ ( ( r > 0 ∧ ( r ( 2 r − 1 ) ( crx 2 − 2 r ) < 0 ∨ ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) > 0 ∧ 6 r < 1 ) ) ) ∨ crx 2 + 4 r 2 = 0 ∨ r 2 ( 2 r − 1 ) ( 16 r 4 − crx 4 ) ≤ 0 ∨ crx 2 r 2 + 4 r 4 > 0 ∨ crx 8 r 2 + 64 r 8 ≤ 4 crx 4 ( 4 r 6 + r 4 ) ) ) ) ∧ ( r = 0 ∨ crx 2 ≥ 4 r 2 ∨ ( ( crx = 0 ∨ crx 2 ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≥ 0 ∨ ( ( crx 2 r ( 6 r − 1 ) ≤ 0 ∨ ( r > 0 ∧ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ≤ 0 ) ) ∧ ( crx 2 r ( 6 r − 1 ) < 0 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 = 0 ) ) ∨ crx 2 r ( 2 r − 1 ) ≥ 0 ∨ ( crx 8 r 2 + 16 crx 2 r 6 ≤ 4 crx 4 ( crx 2 + 1 ) r 4 ∧ ( crx 2 + 4 r 2 = 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) = 0 ) ∧ ( crx 2 ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ( crx 2 + 4 r 2 ) 2 ≥ 0 ∨ ( crx = 2 r ∧ crx 2 = 2 r ∧ crx + 2 r = 0 ∧ crx 2 + 2 r = 0 ) ) ) ) ∧ ( 6 r = 1 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ≥ 0 ∨ r ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) ) > 0 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 8 r 2 ( 6 r ( 3 r − 1 ) + 1 ) < 0 ) ∧ ( crx = 2 r ∨ ( 0 < r < 6 1 ∧ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 8 r 2 ( 6 r ( 3 r − 1 ) + 1 ) > 0 ) ∨ crx + 2 r = 0 ∨ crx 2 + 4 r 2 = 0 ∨ ( crx 2 − 1 ) r 2 ( crx 4 − 16 r 4 ) ≥ 0 ) ∧ ( crx = 2 r ∨ 0 < r < 6 1 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 8 r 2 ( 6 r ( 3 r − 1 ) + 1 ) < 0 ∨ crx + 2 r = 0 ∨ crx 2 + 4 r 2 = 0 ∨ ( crx 2 − 1 ) r 2 ( crx 4 − 16 r 4 ) ≥ 0 ) ∧ ( crx = 2 r ∨ r ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) ) < 0 ∨ crx + 2 r = 0 ∨ crx 2 + 4 r 2 = 0 ∨ 0 < r < 6 1 ∨ r 2 ( 2 r − 1 ) ( 16 r 4 − crx 4 ) < 0 ∨ 4 crx 4 ( 4 r 6 + r 4 ) ≤ crx 8 r 2 + 64 r 8 ) ∧ ( crx = 2 r ∨ r ( 2 r − 1 ) ( crx 2 − 2 r ) < 0 ∨ crx + 2 r = 0 ∨ crx 2 + 4 r 2 = 0 ∨ ( crx 2 ( 3 − 10 r ) < 2 r ( 2 r ( 8 r − 3 ) + 1 ) ∧ 6 r < 1 ) ∨ r 2 ( 2 r − 1 ) ( 16 r 4 − crx 4 ) < 0 ∨ 4 crx 4 ( 4 r 6 + r 4 ) ≤ crx 8 r 2 + 64 r 8 ) ∧ ( 6 r < 1 ∨ 2 r 2 ≥ r ∨ ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≥ 0 ∨ ( ( crx 2 + 4 r 2 = 0 ∨ crx 4 + 4 crx 2 r ( 3 r − 2 ) + 8 r 2 ( 8 r 2 − 4 r + 1 ) = 0 ) ∧ ( ( crx = 2 r ∧ crx + 2 r = 0 ) ∨ crx 2 r 2 + 4 r 4 ≤ 0 ) ) ∨ ( r 2 ( 2 r − 1 ) ( 4 r 2 − crx 2 ) ≤ 0 ∧ ( crx 2 + 4 r 2 = 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) = 0 ) ) ) ) ) ( 22 )
若在命题 (21) 成立的条件下, 进一步让 (18) 加入大前提, 则期待 (16) 能保证 (11). 用 MMA 验证,
Resolve[ ForAll[ { cx, sx, cy} , cx^ 2 + sx^ 2 == 1 & & - 1 <= cy <= 1 , ( crx* sx) ^ 2 + ( 2 * cx* r + 2 * cy* r - 2 * r - 2 * r + 1 ) ^ 2 <= 1 ] , Reals]
FullSimplify[ %, 0 < Abs [ crx] <= 2 r & & Abs [ - 2 * r - 2 * r + 1 ] + 2 r + 2 r <= 1 ]
得
Cr x 2 + 36 r 2 ≤ 12 r . (23) \tag{23}
\text{Cr}_x^2 + 36 r^2 \leq 12 r. Cr x 2 + 36 r 2 ≤ 12 r . ( 23 )
即 MMA 认为, 在 (16)(18)(21) 成立的大前提下, 命题 (11) 等价于命题 (23), 而后者并非恒成立; 这就否定了 (16) 对 (11) 的充分性, 表明 MMA 不支持本节的主要结果 (16). 那么, 问题出在哪? 请求读者帮助! 文末附有联系方式 .
Problem 1.3
本节讨论格式 (3) 的精度. 格式 (3a) 的 modified equation
B t + u ⋅ ∇ B − ∇ ⋅ ( D ∇ B ) + K r B = − τ 2 B t 2 − h x 2 6 u B x 3 − h y 2 6 v B y 3 + h x 2 12 D x B x 4 + h y 2 12 D y B y 4 + O ( τ 2 + h x 2 + h y 2 ) . \begin{equation*}
\tag{24}
\begin{aligned}
& \qquad B_t + \boldsymbol{u} \cdot \nabla B - \nabla \cdot (D \nabla B) + K_{\text{r}} B \\
& = -\frac{\tau}{2} B_{t^2} - \frac{h_x^2}{6} u B_{x^3} - \frac{h_y^2}{6} v B_{y^3} + \frac{h_x^2}{12} D_x B_{x^4} + \frac{h_y^2}{12} D_y B_{y^4} + \mathcal{O}(\tau^2 + h_x^2 + h_y^2).
\end{aligned}
\end{equation*} B t + u ⋅ ∇ B − ∇ ⋅ ( D ∇ B ) + K r B = − 2 τ B t 2 − 6 h x 2 u B x 3 − 6 h y 2 v B y 3 + 12 h x 2 D x B x 4 + 12 h y 2 D y B y 4 + O ( τ 2 + h x 2 + h y 2 ) . ( 24 )
试图把 (24) 右端的时间误差项主部系数 B t 2 B_{t^2} B t 2 用空间导数表示, 所以让 (24) 分别被 ∂ t \partial_t ∂ t 和 ∇ \nabla ∇ 算子作用, 得
B t 2 + u ⋅ ∇ B t − ∇ ⋅ ( D ∇ B t ) + K r B t = O ( τ + h x 2 + h y 2 ) , \begin{equation*}
\tag{25}
B_{t^2} + \boldsymbol{u} \cdot \nabla B_t - \nabla \cdot (D \nabla B_t) + K_{\mathrm{r}} B_t = \mathcal{O}(\tau + h_x^2 + h_y^2),
\end{equation*} B t 2 + u ⋅ ∇ B t − ∇ ⋅ ( D ∇ B t ) + K r B t = O ( τ + h x 2 + h y 2 ) , ( 25 )
∇ B t + ( u ⋅ ∇ ) ∇ B − ∇ [ ∇ ⋅ ( D ∇ B ) ] + K r ∇ B = O ( τ + h x 2 + h y 2 ) . \begin{equation*}
\tag{26}
\nabla B_{t} + (\boldsymbol{u} \cdot \nabla) \nabla B - \nabla \left[ \nabla \cdot (D \nabla B) \right] + K_{\mathrm{r}} \nabla B = \mathcal{O}(\tau + h_x^2 + h_y^2).
\end{equation*} ∇ B t + ( u ⋅ ∇ ) ∇ B − ∇ [ ∇ ⋅ ( D ∇ B ) ] + K r ∇ B = O ( τ + h x 2 + h y 2 ) . ( 26 )
由 (25)(26) 和 (24) 得
B t 2 = ( u ⋅ ∇ ) [ ( ( u ⋅ ∇ ) − K r ) B − ∇ ⋅ ( D ∇ B ) ] − K r [ u ⋅ ∇ B − ∇ ⋅ ( D ∇ B ) + K r B ] − ∇ ⋅ [ D ( u ⋅ ∇ ) ∇ B − D ∇ [ ∇ ⋅ ( D ∇ B ) ] + K r D ∇ B ] + O ( τ + h x 2 + h y 2 ) . \begin{equation*}
\tag{27}
\begin{aligned}
B_{t^2} & = (\boldsymbol{u} \cdot \nabla) \left[ \left( (\boldsymbol{u} \cdot \nabla) - K_{\text{r}} \right) B - \nabla \cdot (D \nabla B) \right] \\
& - K_{\mathrm{r}} \left[ \boldsymbol{u} \cdot \nabla B - \nabla \cdot (D \nabla B) + K_{\text{r}} B \right] \\
& - \nabla \cdot \left[ D (\boldsymbol{u} \cdot \nabla) \nabla B - D \nabla \left[ \nabla \cdot (D \nabla B) \right] + K_{\mathrm{r}} D \nabla B \right] \\
& + \mathcal{O}(\tau + h_x^2 + h_y^2).
\end{aligned}
\end{equation*} B t 2 = ( u ⋅ ∇ ) [ ( ( u ⋅ ∇ ) − K r ) B − ∇ ⋅ ( D ∇ B ) ] − K r [ u ⋅ ∇ B − ∇ ⋅ ( D ∇ B ) + K r B ] − ∇ ⋅ [ D ( u ⋅ ∇ ) ∇ B − D ∇ [ ∇ ⋅ ( D ∇ B ) ] + K r D ∇ B ] + O ( τ + h x 2 + h y 2 ) . ( 27 )
Eq. (27) 有些复杂, 感兴趣的读者可尝试将其代入 (24), 整理出只含空间导数的截断误差项主部, 看有无可能通过适当选取网格参数 τ , h x , h y \tau, h_x, h_y τ , h x , h y 来消去截断误差中的 O ( τ ) \mathcal{O}(\tau) O ( τ ) 项. 我对此不抱多少期望, 所以我说格式 (3a) 是 O ( τ + h x 2 + h y 2 ) \mathcal{O}(\tau + h_x^2 + h_y^2) O ( τ + h x 2 + h y 2 ) 精度的, so does (3b), and their combination (3).
Problem 1.4
将格式 (3) 写成迭代形式,
[ B O ] i , j n + 1 = ( 1 − 2 r x − 2 r y ) [ B O ] i , j n − τ [ K r K r K a ] [ B O ] i , j n + τ [ q i , j K a O s a t ] + ∑ k ∈ { − 1 , 1 } [ ( r x − C r x 2 s g n k ) [ B O ] i + k , j n + ( r y − C r y 2 s g n k ) [ B O ] i , j + k n ] , \begin{equation*}
\tag{28}
\begin{aligned}
\begin{bmatrix*}
B \\
O
\end{bmatrix*}^{n+1}_{i,j} & = (1 - 2 r_x - 2 r_y) \begin{bmatrix*}
B\\
O
\end{bmatrix*}^n_{i,j} - \tau \begin{bmatrix*}
K_{\mathrm{r}} & \\
K_{\mathrm{r}} & K_{\mathrm{a}}
\end{bmatrix*} \begin{bmatrix*}
B \\
O
\end{bmatrix*}^{n}_{i,j} + \tau \begin{bmatrix}
q_{i,j} \\
K_{\mathrm{a}} O_{\mathrm{sat}}
\end{bmatrix} \\
& + \sum_{k \in \{ -1, 1 \}} \left[
\left( r_x - \frac{\mathrm{Cr}_x}{2} \mathop{\mathrm{sgn}}{k} \right) \begin{bmatrix}
B \\
O
\end{bmatrix}^n_{i+k, j} + \left( r_y - \frac{\mathrm{Cr}_y}{2} \mathop{\mathrm{sgn}}{k} \right) \begin{bmatrix}
B \\
O
\end{bmatrix}^n_{i, j+k}
\right],
\end{aligned}
\end{equation*} [ B O ] i , j n + 1 = ( 1 − 2 r x − 2 r y ) [ B O ] i , j n − τ [ K r K r K a ] [ B O ] i , j n + τ [ q i , j K a O sat ] + k ∈ { − 1 , 1 } ∑ [ ( r x − 2 Cr x sgn k ) [ B O ] i + k , j n + ( r y − 2 Cr y sgn k ) [ B O ] i , j + k n ] , ( 28 )
其中 BOD 的源项
q i , j = ∑ ( x 0 , y 0 ) Q B ( x 0 , y 0 ) h x h y I i , j ( x 0 , y 0 ) , \begin{equation*}
\tag{29}
q_{i,j} = \sum_{(x_0, y_0)} \frac{Q_{\mathrm{B}}(x_0, y_0)}{h_x h_y} \mathrm{I}_{i,j} (x_0, y_0),
\end{equation*} q i , j = ( x 0 , y 0 ) ∑ h x h y Q B ( x 0 , y 0 ) I i , j ( x 0 , y 0 ) , ( 29 )
where ( x 0 , y 0 ) (x_0, y_0) ( x 0 , y 0 ) 是 BOD 的一个源点, Q B ( x 0 , y 0 ) Q_{\mathrm{B}}(x_0, y_0) Q B ( x 0 , y 0 ) 是这源的强度 (加载率), I i , j ( x 0 , y 0 ) \mathrm{I}_{i, j}(x_0, y_0) I i , j ( x 0 , y 0 ) 是依 (4) 定义的示性函数. 在网格边界, 按 (7)(8') 设置 halo points. 用 Matlab 实现数值求解二维 Streeter-Phelps 模型 (1)(2) 的 FTCS 方案 (28).
在 (28) 中, 令
[ B O ] n + 1 = [ B O ] n = : [ B ∗ O ∗ ] , \begin{equation*}
\tag{30}
\begin{bmatrix*}
B \\
O
\end{bmatrix*}^{n+1} = \begin{bmatrix*}
B \\
O
\end{bmatrix*}^{n} =: \begin{bmatrix*}
B^{*} \\
O^{*}
\end{bmatrix*},
\end{equation*} [ B O ] n + 1 = [ B O ] n =: [ B ∗ O ∗ ] , ( 30 )
得数值稳态 (steady-state) 解 B ∗ , O ∗ B^*, O^* B ∗ , O ∗ 满足的方程
τ [ q i , j K a O s a t ] = ( 2 r x + 2 r y ) [ B ∗ O ∗ ] i , j + τ [ K r K r K a ] [ B ∗ O ∗ ] i , j − ∑ k ∈ { − 1 , 1 } [ ( r x − C r x 2 s g n k ) [ B ∗ O ∗ ] i + k , j + ( r y − C r y 2 s g n k ) [ B ∗ O ∗ ] i , j + k ] . \begin{equation*}
\tag{31}
\begin{aligned}
\tau \begin{bmatrix}
q_{i,j} \\
K_{\mathrm{a}} O_{\mathrm{sat}}
\end{bmatrix} & = (2 r_x + 2 r_y) \begin{bmatrix*}
B^{*} \\
O^{*}
\end{bmatrix*}_{i,j} + \tau \begin{bmatrix*}
K_{\mathrm{r}} & \\
K_{\mathrm{r}} & K_{\mathrm{a}}
\end{bmatrix*} \begin{bmatrix*}
B^{*} \\
O^{*}
\end{bmatrix*}_{i,j} \\
& - \sum_{k \in \{ -1, 1 \}} \left[
\left( r_x - \frac{\mathrm{Cr}_x}{2} \mathop{\mathrm{sgn}}{k} \right) \begin{bmatrix}
B^{*} \\
O^{*}
\end{bmatrix}_{i+k, j} + \left( r_y - \frac{\mathrm{Cr}_y}{2} \mathop{\mathrm{sgn}}{k} \right) \begin{bmatrix}
B^{*} \\
O^{*}
\end{bmatrix}_{i, j+k}
\right].
\end{aligned}
\end{equation*} τ [ q i , j K a O sat ] = ( 2 r x + 2 r y ) [ B ∗ O ∗ ] i , j + τ [ K r K r K a ] [ B ∗ O ∗ ] i , j − k ∈ { − 1 , 1 } ∑ [ ( r x − 2 Cr x sgn k ) [ B ∗ O ∗ ] i + k , j + ( r y − 2 Cr y sgn k ) [ B ∗ O ∗ ] i , j + k ] . ( 31 )
要求解 (31), 可先求解 (31) 中关于 { B i , j ∗ } \{ B^*_{i,j} \} { B i , j ∗ } 的常系数线性方程. 引入 halo points (7)(8') 后, 这方程的系数矩阵是确定的方阵, 且含大量零元素 (是 sparse matrix (稀疏矩阵)), 适合利用 Matlab 的 sparse() 函数来存储. 求出 B ∗ B^* B ∗ 后, 从 (31) 和 (7)(8') 可建立关于 { O i , j ∗ } \{ O^*_{i,j} \} { O i , j ∗ } 的常系数线性方程, 从而确定 FTCS 格式 (28) 的数值稳态解 (30). 我们期待数值解有渐近行为
[ B O ] n → [ B ∗ O ∗ ] , ( n → ∞ ) . \begin{equation*}
\tag{32}
\begin{bmatrix*}
B \\
O
\end{bmatrix*}^n \to \begin{bmatrix*}
B^{*} \\
O^{*}
\end{bmatrix*}, \quad (n \to \infty).
\end{equation*} [ B O ] n → [ B ∗ O ∗ ] , ( n → ∞ ) . ( 32 )
为确保 FTCS 格式 (28) 稳定且收敛, 应适当选取网格参数 h x , h y h_x, h_y h x , h y 和 τ \tau τ . 根据第 1.3 节的分析, 当 v = 0 , D x = D y = D v = 0, \, D_x = D_y = D v = 0 , D x = D y = D , 格式稳定的一个充分不必要条件 (19) 等价于
h x = h y = : h ≤ 2 D / ∣ u ∣ , τ ≤ ( 4 D / h 2 + max { K a , K r } ) − 1 . \begin{equation*}
\tag{33}
h_x = h_y =: h \le {2 D / |u|}, \quad \tau \le \left(4 D / h^2 + \max{\{ K_{\mathrm{a}}, K_{\mathrm{r}} \}} \right)^{-1}.
\end{equation*} h x = h y =: h ≤ 2 D /∣ u ∣ , τ ≤ ( 4 D / h 2 + max { K a , K r } ) − 1 . ( 33 )
由 (33), 我们选择 h = 5 × 1 0 − 2 m , τ = 1.24 × 1 0 − 3 s h = 5 \times 10^{-2} \, \text{m}, \, \tau = 1.24 \times 10^{-3} \, \text{s} h = 5 × 1 0 − 2 m , τ = 1.24 × 1 0 − 3 s .
我们测试了以下四种排放方案:
单点排放 (single-point discharge). 所有 BOD 都排放在 ( x , y ) = ( 10 , 15 ) m (x, y) = (10, 15) \, \text{m} ( x , y ) = ( 10 , 15 ) m 处.
两点排放 (two-points discharge). BOD 被平均地排放在 ( x , y ) = ( 10 , 15 ) , ( 70 , 15 ) (x, y) = (10, 15), (70, 15) ( x , y ) = ( 10 , 15 ) , ( 70 , 15 ) 两个点上.
线分布排放 (line-distributed discharge). BOD 排放源均匀分布在线段 x ∈ [ 10 , 70 ] , y = 15 x \in [10, 70], \, y = 15 x ∈ [ 10 , 70 ] , y = 15 上.
面分布排放 (face-distributed discharge). BOD 排放源均匀分布在区域 ( x , y ) ∈ [ 10 , 70 ] × [ 10 , 20 ] (x, y) \in [10, 70] \times [10, 20] ( x , y ) ∈ [ 10 , 70 ] × [ 10 , 20 ] 上.
测试结果示于下面 8 个图. 每个图都标注了相应时刻的 BOD 或 DO 的最值. 发现以下特点:
对于点排放 (单点, 两点) 情形, BOD 最高值逐渐升高至稳态值, 且位于最下游的排放源点. DO 最低值逐渐下降至稳态值, 且位置随时间推移有向下游移动的趋势, 最终稳定在排放源点的下游.
对于均匀分布排放 (线, 面) 情形, BOD 最高值逐渐升高至稳态值, 且位置随时间推移有从源区中心向下游移动的趋势, 最终稳定在源区下游边界的内侧. DO 最低值逐渐下降至稳态值. 随时间推移, DO 最低值点有从源区中心向下游移动的趋势, 且相对于 BOD 最高值点, 有从滞后 (即位于相对上游) 到超前的移动趋势.
采取分布排放, 可以大幅降低 BOD 峰值, 因为相比点排放可以大幅降低排放源强度. 但未必能使 DO 最低值达到最高. 物理解释: 点排放引起 DO 在小范围内有快速下降的倾向, 分布排放引起 DO 在大范围内有缓慢下降的倾向. 若平流速度充分大或点排放源强度相对足够小, 则对于点排放情形, 水团受点排放源的影响时间较小, 引起的 DO 损失足够被水-气氧气通量来补偿; 而对于分布排放, 水团在较长时间内持续地被 BOD 反应掉, 足以产生比点排放情形更低的 DO 最低值.